home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Info-Mac 4
/
Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso
/
Science
/
RLaB
/
toolbox
/
rk4.r
< prev
next >
Wrap
Text File
|
1994-04-25
|
1KB
|
78 lines
//
// RLaB 4th Order Runge-Kutta Integrator
// Fixed step size.
// Ian Searle, 6/22/92
//
// Inputs to rk4()
// eval: User-Function to evaluate state-derivative.
// has arguments (time, y)
// tzero: Start time
// tfinal: Finish time
// y0: Initial conditions
// dt: Integration time step
// Output: Matrix of time, and state vector.
//
static (collapse, rrk4);
rk4 = function(eval, tzero, tfinal, y0, dt)
{
local(i, y, lout, I);
y = y0[:];
lout = << [tzero, y0.'] >>; I = 2;
for(i in tzero:(tfinal-dt):dt)
{
lout.[I] = [i+dt, (y = rrk4 (y, i, dt, eval))'];
I++;
}
return collapse (lout);
};
rrk4 = function(y, x, h, eval)
{
local(i, xh, hh, h6, dydx, dym, dyt, yt, yout);
// Initialize
hh = h*0.5;
h6 = h/6;
xh = x + hh;
dydx = eval(x, y); // 1st step
yt = y + hh*dydx;
dyt = eval(xh, yt); // 2nd step
yt = y + hh*dyt;
dym = eval(xh, yt); // 3rd step
yt = y + h*dym;
dym = dym + dyt;
dyt = eval(x+h, yt); // 4th step
yout = y + h6*(dydx + dyt + 2.0*dym);
return yout;
};
//
// Collapse a LIST of matrices into a single matrix.
//
collapse = function ( LIST )
{
local (i, m, row, col);
row = size (LIST.[members (LIST)[1]])[1];
col = size (LIST.[members (LIST)[1]])[2];
m = zeros (length (LIST)*row, col);
for (i in 1:length (LIST))
{
m[i;] = LIST.[i];
}
return m;
};